# ===============================================
# Paris ratchet paper
# Code for summary statistics and plots
# ===============================================


## ## Contents

#' 1. Load in data for models from "00_make_analysis_data.R"
#' 2. Summary statistics
#' 3. UK emissions and targets plot
#' 4. Paris v. Glasgow scatterplot
#' 5. Ratchet map
#' 6. Descriptive statistics on ratchet


## Objects

#### Initialize #### 

library(tidyverse)
library(knitr)
library(DescTools) # Winsorize()
library(broom)
library(modelsummary)
library(nnet) # for multinomial logit
library(magrittr)
library(countrycode)


# Set working directory

country_level_average <- readRDS("output/country_level_average_20230920.rds")
data_for_models <- readRDS("output/data_for_models20230920.Rds")


#### How many ratchet? ####

#' These are the facts in-text in section 5

## What share of global emissions in ratchet
{country_level_average %>% 
  filter(year == 2019) %>% 
  filter(!is.na(ndc2_date)) %>% 
  summarize(sum(ghg_cw_ex))} /
  {country_level_average %>% 
      filter(year == 2019) %>% 
      # filter(!is.na(ndc2_date)) %>% 
      summarize(sum(ghg_cw_ex))} 
#' 94% of emissions

# Directions
data_for_models %>% 
  mutate(ratchet = case_when(glasgow_target_safe > paris_target_safe ~ "Upgrade",
                             glasgow_target_safe == paris_target_safe ~ "Constant",
                             glasgow_target_safe < paris_target_safe ~ "Weaken")) %>% 
  group_by(ratchet) %>% 
  summarize(n())

# Average change within groups
data_for_models %>% 
  mutate(ratchet = case_when(glasgow_target_safe > paris_target_safe ~ "Upgrade",
                             glasgow_target_safe == paris_target_safe ~ "Constant",
                             glasgow_target_safe < paris_target_safe ~ "Weaken")) %>% 
  group_by(ratchet) %>% 
  summarize(median_diff = median(glasgow_target_safe - paris_target_safe))

### Summary statistics

table_summary_statistics <- data_for_models %>% 
  select(paris_target_safe, glasgow_target_safe,
         io_percentage,
         tradeflow_percentage, eite_percentage, 
         competition_tradeflow_percentage, competition_eite_percentage) %>% 
  filter(!is.na(paris_target_safe),
  !is.na(glasgow_target_safe)) %>%
  pivot_longer(names_to = "variable", values_to = "value",
               paris_target_safe:competition_eite_percentage) %>% 
  filter(!is.na(value)) %>% 
  group_by(variable) %>% 
  summarize(observations = n(),
            # missing = 192 - observations,
            mean = mean(value),
            sd = sd(value),
            min = min(value),
            median = median(value),
            max = max(value)) %>% 
  mutate_at(vars(mean:max), ~round(., 1)) %>% 
  kable(., format = 'latex', booktabs = T)
table_summary_statistics


# Highly correlated variables
table_correlations <- data_for_models %>% 
  select(glasgow = glasgow_target_safe,
         paris = paris_target_safe,
         io = io_percentage,
         total_trade = tradeflow_percentage,
         eite_trade = eite_percentage) %>% 
  drop_na() %>%
  cor() %>% 
  as.data.frame() %>% 
  mutate_all(~round(., digits = 3)) %>% 
  kable(., format = 'latex', booktabs = T)
table_correlations

#### UK emissions

gbr_to_plot <- country_level_average %>% 
  filter(iso3c == "GBR") %>% 
  select(iso3c, year, ghg_cw_chosen, ndc1_absolute_average, ndc2_absolute_average) %>% 
  add_row(iso3c = "GBR", year = 2030, ghg_cw_chosen = 417.7) %>% 
  add_row(iso3c = "GBR", year = 2030, ghg_cw_chosen = 234.8) %>% 
  add_row(iso3c = "GBR", year = 2050, ghg_cw_chosen = 0) %>% 
  select(iso3c, year, ghg_cw_chosen) %>% 
  mutate(series = c(rep("Observed", times = 30), "NDC 2015", "NDC 2021", "Net zero"),
         series2 = c(rep("Observed", times = 30), "Paris target", "Glasgow target", "Net zero")) 

gbr_ratchet <- ggplot() +
  geom_line(data = gbr_to_plot %>% 
              filter(year < 2020),
            aes(year, ghg_cw_chosen),
            size = 1.5) +
  geom_point(data = gbr_to_plot %>%
              filter(year > 2018),
            aes(year, ghg_cw_chosen)) +
  geom_text(data = gbr_to_plot %>% 
              filter(year > 2020),
            aes(year, ghg_cw_chosen, label = series2, family = "serif"), 
            size = 5,
            hjust = -0.1) +
  geom_line(data = gbr_to_plot %>% 
              filter(series %in% c("Observed", "NDC 2015"),
                     year > 2018),
            aes(year, ghg_cw_chosen),
            linetype = 2) +
  geom_line(data = gbr_to_plot %>% 
              filter(series %in% c("Observed", "NDC 2021"),
                     year > 2018),
            aes(year, ghg_cw_chosen),
            linetype = 2) +
  geom_line(data = gbr_to_plot %>% 
              filter(series %in% c("Observed", "Net zero"),
                     year > 2018),
            aes(year, ghg_cw_chosen),
            linetype = 2) +
  scale_y_continuous(limits = c(0, 750),
                     expand = c(0.1,0),
                     breaks = seq(0, 750, 150)) +
  scale_x_continuous(limits = c(1990, 2060),
                     breaks = seq(1990, 2050, 10)) +
  labs(x = "", y = "Annual GHG (MtCO2)") +
  theme(legend.position = "bottom", 
        panel.background = element_rect(fill=NA),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y =  element_blank(),
        panel.grid.minor.x =  element_blank(),
        panel.ontop = FALSE,
        axis.line = element_line(size = 1, colour = "black"),
        axis.text = element_text(size = 12, family = "serif"),
        strip.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        axis.title = element_text(size = 12, family = "serif"))
ggsave(plot = gbr_ratchet, 
       filename = "plots/gbr_ratchet.pdf", 
       units = "in", height = 4, width = 6)


#### Ratchet scatterplot #### 

paris_v_glasgow <- data_for_models %>% 
  select(iso3c, paris_target_safe, glasgow_target_safe) %>% 
  mutate(gbr_flag = ifelse(iso3c == "GBR", T, F)) %>% 
  ggplot(., aes(paris_target_safe, glasgow_target_safe, 
                label = iso3c)) +
  scale_x_continuous(breaks = seq(-200, 100, 50)) +
  scale_y_continuous(breaks = seq(-200, 100, 50)) +
  coord_cartesian(xlim = c(-215, 100), ylim = c(-215, 100)) +
  geom_abline(slope = 1, intercept = 0, linetype = 1, color = 'grey80') +
  # geom_hline(yintercept = 0, linetype = 2, color = 'grey90') +
  # geom_vline(xintercept = 0, linetype = 2, color = 'grey90') +
  geom_point(aes(color = gbr_flag, shape = gbr_flag, size = gbr_flag)) +
  scale_shape_manual(values = c(1, 19)) +
  scale_color_manual(values = c("black", "#1b7837")) +
  scale_size_manual(values = c(2, 2.5)) +
  labs(color = NULL, shape = NULL, size = NULL,
       x = "Paris GHG target", y = "Glasgow GHG target") +
  theme(legend.position = "none", 
        panel.background = element_rect(fill=NA),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y =  element_blank(),
        panel.grid.minor.x =  element_blank(),
        axis.line = element_line(size = 1, colour = "black"),
        axis.text = element_text(size = 12, family = "serif"),
        strip.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        axis.title = element_text(size = 12, family = "serif"),
        panel.ontop = FALSE)
ggsave(plot = paris_v_glasgow, 
       filename = "plots/paris_v_glasgow.pdf", 
       units = 'in', height = 4, width = 6)


#### Map ####


WorldData <- map_data('world')
WorldData <- WorldData 
WorldData <- fortify(WorldData)
WorldData$iso3c <- countrycode(WorldData$region, origin = 'country.name', destination = 'iso3c')


ratchet_map_df <- data_for_models %>% 
  select(iso3c, paris_target_safe, glasgow_target_safe, paris_target_impute) %>% 
  mutate(ratchet_amount = glasgow_target_safe - paris_target_safe,
         ratchet_binary = ifelse(glasgow_target_safe > paris_target_safe, 1, 0),
         ratchet_bucket = case_when(ratchet_amount < -50 ~ "Weaker than -50",
                                    ratchet_amount >= -50 & ratchet_amount < -25 ~ "-50 to -25",
                                    ratchet_amount >= -25 & ratchet_amount < 0 ~ "-25 to -1",
                                    ratchet_amount == 0 ~ "0",
                                    ratchet_amount > 0 & ratchet_amount < 25 ~ "1 to 25",
                                    ratchet_amount >= 25 ~ "Over 25")) %>% 
  mutate(ratchet_bucket = factor(ratchet_bucket),
         ratchet_bucket = ordered(ratchet_bucket, levels = c("Weaker than -50", 
                                                             "-50 to -25",
                                                             "-25 to -1",
                                                             "0",
                                                             "1 to 25",
                                                             "Over 25"))) %>% 
  full_join(., WorldData, by = "iso3c") %>% 
  arrange(order)

## Purple to green color scale
my_fill_scale <- c('#762a83','#9970ab','#c2a5cf','grey80', "#a6dba0", "#1b7837")

# Discrete color scale
ratchet_map <- ggplot() +
  geom_map(data = WorldData, 
           map = WorldData,
           aes(x = long, y = lat, 
               group = group, map_id = region),
           fill = "black", 
           colour = "#7f7f7f", size = 0) + 
  geom_map(data = ratchet_map_df, 
           map = ratchet_map_df,
           aes(fill = ratchet_bucket, map_id = region),
           colour ="#7f7f7f", size = 0) +
  coord_map("rectangular", lat0 = 0, 
            xlim = c(-180,180), 
            ylim = c(-60, 90)) + 
  scale_fill_manual(values = my_fill_scale, na.value = "white") +
  scale_y_continuous(breaks = c(), expand = c(0,0)) +
  scale_x_continuous(breaks = c(), expand = c(0,0)) +
  labs(fill = "Ratchet percentage", 
       title=NULL, x=NULL, y=NULL) +
  theme_void() +
  guides(fill = guide_legend(nrow = 2, byrow = T)) +
  theme(panel.border = element_blank(),
        text = element_text(size=12, family = "serif"),
        legend.key.size = unit(0.25, 'cm'),
        legend.position = "bottom") #
ggsave(plot = ratchet_map, 
       filename = "plots/ratchet_map_discrete.pdf", units = "in", height = 5, width = 9)










